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Abstract. We present theoretical and experimental results on high-fidelity 
transfer of a trapped Bose-Einstein condensate into its first vibrationally excited 
eigenstate. The excitation is driven by mechanical motion of the trap, along a 
trajectory obtained from optimal control theory. Excellent agreement between 
theory and experiment is found over a large range of parameters. We develop 
an approximate model to map the dynamics of the many-body condensate wave 
function to a driven two- level system. 
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1. Motivation 

The precise control over quantum systems represents a major challenge in modern 
physics. Successful implementation of quantum technologies may lead to the 
construction of devices such as quantum simulators, quantum cryptography devices, 
and quantum computers. For such applications, one needs to produce arbitrary 
quantum states, e.g. strongly entangled many-body states, or states which are far 
from thermal equilibrium or the ground state of the system. 

In this article we report on highly efficient preparation of a non-classically 
excited motional state of a Bose-Einstein condensate (BEC), by a modulation of the 
trapping potential, as obtained from optimal control theory (OCT). Fast changes 
of the potential are routinely used in BEC laboratories, for instance as ways to 
probe the gas by exciting collective excitations [HE], or to displace the samples 
for further manipulation. Recently, controlled modulations of the trapping potential 
were achieved in order to quickly displace BECs while keeping them in their ground 
state. These "shortcuts to adiabaticity" [31 [4| [5] |6j |7] take advantage of the self-similar 
dynamics of interacting BECs trapped in time-dependent harmonic potentials [8^, '9] . 

For more general desired states, like excited stationary states, for which no exact 
solutions are found, one needs to use numerical methods such as OCT. Such approaches 
were investigated theoretically for the splitting of BECs [lOj, to optimize the transport 
of atoms in optical lattices for quantum gate operations [llj, or to amplify number 
squeezing [12j. 

Here, we aim for a vibrational state inversion, where the entire population of the 
condensate is transferred to the first excited state of its motional degree of freedom. 
Such an inverted state can be used as a source for the amplified emission of matter- 
wave twin beams [TF, 14J , similar to a pumped gain medium in a laser or an optical 
parametric amplifier [Fig.[TJa)].f We start from a condensate in the ground state along 
the strongly confined (transverse) directions of an elongated trapping potential. We 
then use OCT on a controlled displacement of the trap center [transverse "shaking" 
of the cloud. Fig. [ijc)], in order to transfer the BEC to the first antisymmetric 
stationary state as given by the Gross-Pitaevskii equation (GPE), which is governing 
the system's dynamics [Fig.[TJb)]. The efficiency of this process is close to 100%, which 
corresponds to the desired vibrational inversion of the atomic cloud. Since a time- 
dependent harmonic potential (where all energy levels are equidistant) would not allow 
to transfer to an excited stationary state fTS] [161 HZl 'TS] , we here use an anharmonic 
potential [19^.20] generated by a radio-frequency dressed magnetic trap [21j. For a 
non-interacting gas, the final state would simply correspond to all the atoms residing 
in the first excited eigenstate of the trap (quantum numbers = 0, ny = l,nz = 0). 
However, in our many-body wave function inter- atomic interactions are an essential 
ingredient of the system's dynamics and cannot be neglected in the optimization. In 
the experiment they manifest themselves in energy shifts due to the atomic mean-field, 
and a decay of the excited state by means of inelastic two-body scattering [13, 14J. 

Another aspect we will address is the interpretation of our results beyond a 
simple comparison of calculated and measured wave function dynamics. While such a 
comparison benchmarks the accuracy to which experiments and theory are matched, 
it provides only limited insight into the nature of the excitation mechanism, and the 
structure of the quantum state during and after the excitation process. To this end, 
we will deduce an approximate description, that allows to map the many-body wave 
function in a weakly anharmonic confinement to a driven two-level system, where the 
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Figure 1. Schematic representation of our excitation scheme, (a) Ihustration of 
the sequence: an elongated Bose-Einstein condensate (blue) is trapped on an atom 
chip (chip and condensate aspect ratio not drawn to scale). Mechanical motion 
along the y-axis pumps the gas into a vibrationally excited state, which decays 
by directed emission of twin-atom beams along x (red), (b) Initial and excited 
condensate wave functions in the (x,y)-plane. The state ilJd{^-> u) (desired state) is 
the first excited state of the Gross-Pitaevskii equation along y. (c) Trajectory 
X(t) of the trap minimum along y. (d) Time-dependent density n(y, t) = \'il^{y^ t)P 
of the condensate wave function under infiuence of the excitation process, (e) 
Population (simulated) of excited state ipd(y) derived from wave- function overlap 
(black, with markers) and using the two-mode model as introduced in Sec. |4.3| 
(red, solid). 

excitation process corresponds to a 7r-pulse that transfers all population to the excited 
state. 

The paper is structured as follows: section [2] presents the theoretical description 
of the problem and the OCT algorithm used to obtain the excitation trajectory of the 
trap center, section |3] details the experimental implementation, and finally, section [4] 
discusses the results, with an emphasis on how the behavior of key observables can be 
captured by a two-level model. To our knowledge, this excitation sequence represents 
the first successful use of OCT for the preparation of exotic many-body states of 
Bose-Einstein condensates. 

2. Optimal control theory 

Optimal control theory (OCT) is a mathematical tool that allows to determine an 
optimal control sequence for a given control problem \22\ |23] . In the following 
we review the basic ingredients of optimal control theory taking the example of 
the shaking process which brings the condensate from the ground state to its first 
vibrationally excited eigenstate. Our analysis closely follows the presentation given in 
Refs. [HHli. 

As will be discussed in Sec. [3j the condensate can be described by a one- 
dimensional Gross-Pitaevskii equation for the transverse coordinate along which 
the condensate is displaced, according to 

.^d^ = + V,{y,t)+gmy,t)f^ ^{y,t) . (1) 



Here m is the mass of the ^^Rb atoms and ^ is a nonlinearity parameter accounting for 
the repulsive atom- atom interactions [24^ 12j. The anharmonic confinement potential 
Vx{y,t) = V6{y — A(t),0) (see Sec. 3.1) follows a control parameter A(t), in our case 
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the displacement of the potential minimum, and provides the means for exciting the 
condensate. The objective of the control problem can now be formulated as follows. 
Let Ao be the control parameter at the initial time t = 0, and Ai the control parameter 
at the final time t = T oi the control process. Likewise, we denote the initial ground 
state of the GPE with '^/^o(^) and the desired final wave function (in our case the first 
excited state of the GPE in the anharmonic trap) with tpdiv)- OCT then seeks for the 
optimal time variation of X{t) that brings the final wave function as close as possible 
to the desired state ip^^. 

To gauge the success of the excitation process for a given control field A(t), we 
define a cost function 



J{i^{T),X) 



i(V'div(T))r 



Jo 



dt. 



(2) 



The first term of the cost function becomes minimal when the final wave function 
precisely matches the desired wave function, apart from a global (irrelevant) phase. 
The second term favors smooth control fields and is needed to make the OCT problem 
well posed [25j. 7 is a parameter that weights the relative importance of the two 
control objectives of smooth control fields and of wave function matching. As our 
experimental implementation allows fast and precise control of A(t), the parameter 7 
can be set such that the control penalization is always much smaller than the first 
term in Eq. OCT is now seeking for an "optimal control" that minimizes the cost 
function J{tlj{T)^ A), under the condition that the final wave function '^(T) has to be 
obtained from the Gross-Pitaevskii equation of Eq. ([T]) with the initial wave function 
ipo{y)' To turn this constrained minimization problem into an unconstrained one, 
within the OCT framework one introduces a Lagrange function 



L(V^,p,A) = J(V^,A) + 3?e 



i^)dt, 



where the adjoint function p{y, t) acts as a generahzed Lagrange parameter. Here and 
in the following we will, for the sake of brevity, often omit parameters y and t. At 
the minimum of J^ip, A) the Lagrange function has a saddle point, where all three 
derivatives 6L/6tp, 5L/5p and 6L/6X must vanish. Performing the usual functional 
derivatives, we obtain after some variational calculation the following optimality 
system: 



ih 



ih 



dt 
dp 
di 



2m dy-^ 

2m dy^ 
dV\ 



+ Vx + 2g\i;\^]p + gi;^p* 



(3a) 



(3b) 
(3c) 



which has to be solved together with the initial condition V^(0) = tpo^ as well as 
with the constraints on the control field A(0) = Aq and A(T) = Ai. To obtain the 
equation for the adjoint function p, we have performed an integration by parts for the 
term involving the time derivative of prior to working out the functional derivative 
SL/Sij). This procedure gives, in addition to Eq. ([3b]), the terminal condition 



(4) 
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Quite generally, the Lagrange parameter determines the sensitivity of the system with 



respect to the external control. In our case, the dynamic equation (3b) describes the 
propagation of fluctuations around the Gross-Pitaevskii solution and is closely related 
to the usual Bogoluibov-de Gennes equations [26j. 

In most cases of interest it is impossible to guess \{t) such that Eqs. ([3^-c) are 
fulfilled simultaneously, and one has to employ a numerical solution scheme. Suppose 



that \{t) is some guess for a viable control field. We can now solve Eq. (3a) forward in 
time to obtain the final wave function V^(T), which, in turn, allows us to compute the 
adjoint function p{T) from Eq. (|4|. In the ensuing step, the time evolution of p(t) is 



solved backwards in time. Since \{t) is not the optimal control, Eq. (3c) is no longer 
fulfilled. However, the functional derivative 

^ = -7A-5Re{^|^b) (5) 

provides us with a search direction for \{t). Adding a fraction of SL/SX to X{t) leads 
to a control that performs better and brings the final wave function '^(T) closer to 
the desired one. The improved control field is then used in the next iteration. In 
our simulations we typically perform a time discretization of the interval [0,T] and 
use a generic optimization routine, such as the nonlinear conjugate gradient [27\ or 
a quasi-Newton method, together with Eq. ^ for computing the appropriate search 
directions. One shortcoming of Eq. ([5| is that in general 5L/SX does not vanish at 
the boundary points of the time interval, although the control field is fixed to the 
values of Aq and Ai there. To overcome this problem, one rewrites the penalization 
term of the control field (7/2) (A, X)^^ (7/2) (A, A)^^!, where the definition of the 
inner product is {u,v)iii = {u^v)l2 [28j. It is important to realize that this different 
norm does neither affect the value of the cost function nor the Gross-Pitaevskii or 
adjoint equations. However, it does affect the equation for the control field in case of 
a non-optimal A(t), which now satisfies a Poisson equation 

The advantages of Eq. ([6| are that the boundary conditions for X{t) are automatically 
fulfilled and that changes due to large values of the second term on the right-hand 
side are distributed, through the solution of the Poisson equation, over the whole time 
interval. In all our OCT calculations we use Eq. ([6| instead of Eq. ([5|. 

Our OCT implementation relies on a numerical optimization routine and a 
differential equation solver. As for the optimization routine, one can use any generic 
code that, starting from some initial guess for the control field, requires a function 
value (the cost function) together with the derivative of the evaluated function 6L/6X 
to compute a new, improved X{t). When using the H^ norm of Eq. ([6| one must 
ensure that all inner products in the generic code are evaluated as {u^ v)y[i rather 
than {u^v)i^2. In general we observed the best performance for the quasi-Newton 
BFGS optimization [29j, which outperforms the nonlinear conjugate gradient method 
for larger number of iterations in the optimization loop. As for the differential equation 
solver, we usually employ a split operator technique [24] because of its robustness and 
simplicity. 

The OCT optimization is schematically depicted in Fig. |2j One starts with some 
initial guess for the control field. In general, the outcome of the OCT loop does not 
depend critically on the initial X{t) and one can use any reasonable guess, such as 
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Figure 2. Schematics of the OCT optimization loop, which starts with an initial 
guess for the control field \{t) associated with the displacement of the minimum 
of the confinement potential. First, the Gross-Pitaevskii equation with the ® 
initial condensate wave function ipo(y) is solved forwards in time, to obtain (D the 
final wave function ip(y, T) at the terminal time T = 5 ms of the control process, 
which in general deviates significantly from the desired, first excited wave function 
ip(i(y). The density plots in the different panels report the time evolution of the 
square moduli of the different functions. From the knowledge of ip{y, T) and 
i^diy) we can co mput e the ® terminal value of T) via Eq. and solve the 
adjoint equation ( |3b| ) backwards in time ®, to finally come up with a new search 
direction for the optimal control field [Eq. ([gJ] that is used in the next iteration 
of the optimization loop. The solid lines superimposed on X{t) in the panel of 
the adjoint equation depict the search directions. The inset ® shows how the 
cost and derivative for a given control decrease with increasing iterations, until 
© an optimal control is obtained. Here X{t) (magnified by a factor 4) steers the 
system from ipo(y) to the desired wave function at the terminal time of the control 
process. 



in our case some interpolating function between the boundary values of Aq = and 
Ai = 0.1 /im at the terminal time T = 5 ms. Next, (D the condensate wave function 
ilj{y, 0) = i^oiy) is set to the ground state V^o(^) of the anharmonic trap, including the 
nonlinear term of the Gross-Pitaevskii equation [24J, and Eq. (3a) is solved forwards 
in time to obtain d) the terminal wave function ?/^(?/,T). For the initial guess of the 
control, ip{y^T) differs significantly from the desired, first excited state i^diy) of the 
anharmonic trap, which has a node in the middle, as can be also inferred from the 
ensuing time evolution where the trap displacement is held constant. From Eq. Q 
we can compute the d) terminal condition for and ® solve the adjoint equation 
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( [3b| ) backwards in time. Finally, the knowledge of the complete history of ij^iy^t) and 
t) allows us to compute the new search direction through Eq. (|6|, and to pass this 
direction to the optimization routine which will come up with a new, improved A(t), 
which can be used in the next iteration of the optimization loop. 

In the inset (D of Fig. [2] we show how the cost function J(?/;, A) and the derivative 
measure \5L/8\\ evolve with increeasing iterations. Note that the "optimal control" 
corresponds to a minimum of the control landscape, associated with a derivative equal 
to zero, but it is generally not guaranteed that also the cost is small there. However, 
there are indications that under quite broad conditions the OCT loop will come up 
with a \{t) that fulfills the control objective of wave function matching almost perfectly 
|3Q| . In our simulations we typically stop after a given number of iterations or when the 
derivative has become sufficiently small. The resulting \(t) sequence is then called the 
optimal control. As can be seen from the solution of the Gross-Pitaevskii equation on 
the in (D, with this control we closely match the desired wave function at the terminal 
time, with a fidelity of |('?/^d|'0(^))P ~ 1 — 3 • 10~^. Up to a global phase, the wave 
function remains stationary for t > T. 

3. Experimental implementation 

The vibrational state control scheme is realized using an ultra-cold Bose gas trapped 
on an atom chip. The experimental procedure is very close to that described in ^13 | il4 j . 
In brief, a laser-cooled cloud of Rubidium-87 atoms in the |F = 1, mp = — 1) Zeeman 
level is loaded into a strongly elongated atom chip wire trap (SU [32] . Using forced 
evaporative cooling, the gas is brought to a temperature close to quantum degeneracy. 
Then, by means of radio- frequency dressing [33 | [2T ] . the external confinement along 
the two tightly trapped axes ^, ^ is deformed from a harmonic to an anharmonic and 
anisotropic potential (see next section for details). In the anharmonic trap, further 
cooling down to a temperature T well below 50 nK is performed. With TV ~ 800 atoms, 
the gas is now in a quasi-condensate |34l|35] regime, where phase fluctuations prevent 
true condensation of the matter wave along the longitudinal (elongated) direction x. 
However, the energy scales correponding to both temperature {k^T ^ hx 550 Hz) and 
atom interactions (chemical potential /i ^ h x 600 Hz) are well below the trap level 
spacing h x 2 kHz) along the transverse directions y^z. Thus, thermal excitations 
in the transverse degrees of freedom are frozen out, and almost all atoms occupy a 
single transverse mode; along y and the system is hence appropriately described by 
a single condensate wave function. In this paper, we are mostly concerned with the 
dynamics along the transverse direction y, and hence neglect the longitudinal mode 
structure of the quasi-condensate. 

Having prepared the system in this way, we apply the control sequence, while 
monitoring the momentum space distribution of the condensate, as will be described 
in the following sections. 

3.1. Trap preparation 

In a quantum harmonic oscillator, all states that can be addressed by simple 
displacement of the potential are quasi-classical coherent states [11]. This statement 
also holds for a harmonically trapped interacting many-body system, where a quasi- 
classical collective oscillation at the trap frequency fully decouples from more complex 
internal dynamics [TTl [T8] . Hence, transferring the condensate population into an 
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Figure 3. (Main figure) Schematic of the atom chip layout (see Ref. [32 for 
details). The waveguide potential is formed by the current through the trap wire 
along —X and a static bias field By^ adding up to quadrupole field (bent arrows). 
An external offset field along B^, perpendicular to the figure plane, defines the 
Larmor frequency at the trap minimum (loffe-Pritchard field configuration). On 
a separate chip layer, currents in broad wires along y (not shown) provide weak 
longitudinal confinement. The radio frequency dressing currents are applied 
to wires (RF) in parallel to the trapping wire, leading to a RF field along z 
(blue arrows). The resulting anisotropic transverse potential is shown as ellipse 
in the center of the quadrupole. Finally, the modulation of the trap position 
is accomplished by a current in an auxiliary wire (M), leading to a magnetic 
field, aligned at ~ 19° with respect to the z axis (red arrow). (Inset) Field 
configuration for trap position modulation. The transverse trap position is defined 
by cancellation of the chip wire field (brown) and the bias field (green). Adding 
a weak field along z (red) tilts the bias field slightly, leading to a horizontal shift 
of the trap minimum. 

excited, stationary state necessitates an anharmonic potential along the displacement 
direction where the decoupling of collective and internal dynamics breaks down. 
Furthermore, to be robust against excitations in the perpendicular direction 
anisotropy in the transverse plane of the potential is required, causing a detuning 
of trap levels between the directions. 

Initially, the loffe-Pritchard field configuration as created by the chip wires (plus 
external offset fields, see Fig. |3| is rotationally symmetric, and provides harmonic 
trapping along the transverse directions z. For the parameters chosen in our 
experiment, the transverse trap frequency is uq = 4.1kHz in both directions, whereas 
the longitudinal frequency is of the order of 30 Hz. To introduce anharmonicity and 
anisotropy, we apply radio- frequency dressing |36l [2T| [33] . Using two chip wires 
running in parallel to the trapping wire as antennae, the atoms are irradiated by 
a radio- frequency (RF) near field with linear magnetic polarization, which is red- 
detuned by tens of kHz with respect to the atomic Larmor frequency due to the static 
magnetic field at the trap center. The RF field adiabatically mixes the Zeeman levels 
of the F = 1 hyperfine manifold, coupling them to dressed states. This gives rise to an 
energy shift that depends on detuning from the Larmor frequency A(r) and coupling 
strength (Rabi frequency) Q{r). Both quantities are position-dependent, the latter 
because of the changing RF polarization with respect to the local magnetic field that 
modulates the coupling strength. In rotating- wave approximation [21j, the resulting 
potential landscape up to a constant is given by: 

V{r)/h = V^(r)2 + A(r)2. 

The dressing is most effective along the direction perpendicular to the RF polarization; 
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Figure 4. Effects of RF dressing on the transverse trapping potential, (a) 
Potential along the y (displacement) direction, as a function of RF Rabi frequency. 
The detuning is Aq = —55 kHz. At dressing strengths above Q^q ~ 180 kHz, 
splitting of the single potential into a double well occurs, (b) Shift of single- 
particle trap levels vs. dressing strength. Solid and dashed lines correspond to 
perpendicular {y) and parallel {z) directions with respect to the RF polarization, 
respectively. Blue: frequency of harmonic part, as defined in Eq. ([7|. Black, red: 
first and second level spacing of single-particle eigenstates. Inset: initial (grey) 
and dressed (black) potential, each with their first three energy levels. The green 
lines in both panels mark the setting used for the experiments. 



in our case, applying a polarization along the vertical axis z leads to a deformation 
mostly along y. In Fig. Qa), the potential along y is shown as a function of dressing 
strength, expressed as coupling near the trap center. At sufficiently strong coupling, 
splitting of the potential into a double well occurs, which is the typical application of 
RF-dressed potentials [SS", "371 |38l [39^, "40] . However, at lower coupling, this technique 
also allows for the introduction of anharmonicity and anisotropy to a single trap, as 
needed for our scheme. In the experiment, we apply a RF field of ~ 0.84 G peak- 
to-peak amplitude, leading to a coupling VLq = 147 kHz, at a frequency red-detuned 
by Aq = —54 kHz with respect to the Larmor frequency near the trap minimum 
(824 kHz). 

The resulting potential is shown as a green line in Fig. [4]^a). Even though the 
rotating- wave approximation holds well for the used dressing strength [41j, the high 
sensitivity of the excitation protocol to the exact potential shape calls for an exact 
calculation by means of a Floquet analysis [42j. Along two transverse directions the 
result can be approximated by a sixth-order polynomial of the form 



2 / \ 4 / \ 6 



Ve{y^z)/h=^-^(f] +a,(f ) + ( f ) (7) 



2 \ly J J 

X 2 / \ 4 / \ 6 



2 \l 



In this expression, the lengths ly^^ = A^h/ {muy^z)/ (27r) correspond to the characteristic 
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length of the harmonic part. The parameters are given by: 

Vy = 1655 Hz; = 2751 Hz (9) 

ay = 78.2 Hz = -69.6 Hz 

= -0.96 Hz = 9.1Hz 

ly = 265 nm Iz = 206 nm. 

Along the sixth-order term is negligibly small, and the description reduces to a 
Duffing oscillator [43j. 

By solving the Schrodinger equation, the single-particle trap levels of the dressed 
potential can be obtained. The first two level spacings along y and z are shown in 
Fig. [4]^b). For the used parameters (as marked by a green line), the initial degeneracy 
of the level spacings is lifted, and we obtain the excitation energies (zero-point energy 
subtracted) [£^10, £^20, ^01, £^02, £11] /h = [1.84, 3.83, 2.58, 5.21] kHz with Eij denoting 
the i-th and j-th state along y and respectively. The relevant level spacings 
along y are given by ui = 1.84 kHz, = 1.99 kHz, the first level spacing along z 
is Uz = 2.58 kHz. From the corresponding eigenfunction along the z direction, and a 
Thomas-Fermi approximation [44J of the longitudinal profile for N = 800 atoms, we 
can estimate the coupling constanl[f]in Eq. ([T]) by averaging diS g = hx 300 Hz pm [45j. 

In the experiment, characterization of the initial harmonic trap is straightforward, 
using radio- frequency spectroscopy and observation of collective oscillations. On the 
other hand, confirming the (calculated) parameters of the anharmonic dressed trap 
with sufficient accuracy is difficult. Instead, we optimize the experimental control 
parameters (RF field strength and detuning) directly, by comparing the response to 
the control ramp to that determined numerically using those trap parameters. Along 
the longitudinal x-axis, the harmonic trap frequency = 16.3 Hz is determined by 
observation of deliberately excited collective modes of the atom cloud. 



3.2. Control of trap motion 

The transverse movement of the potential is accomplished by applying a time- 
dependent current to an auxiliary wire running parallel to the main trapping wire. 
As shown in the inset of Fig. [3) the additional magnetic field along z causes a slight 
tilt of the homogeneous bias field, which is exactly aligned along y initially. The trap 
minimum position, which is given by the point where the bias field cancels that of the 
trapping wire, is displaced along y. Additionally, the ^/-component of the modulation 
field, which changes the magnitude of the bias field, causes a slight proportional 
movement along z. However, as confirmed by two-dimensional simulations, the 
anisotropy of the transversal potential suppresses any significant influence on the 
excitation along y. From numerical simulations of the field geometry, the movement 
of the trap minimum caused by the current can be calculated as 26nm/mA along y 
and 9nm/mA along z. The geometry of all chip wires and homogeneous offset fields 
involved in trapping and modulation is shown in Fig. [3j 

J Note, that we normalize the wave function to 1, not in Eq. (^. Hence, g incorporates the atom 
number. 
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Figure 5. Effect of filtering due to finite electronics bandwidth. A;-axes 
are scaled to hk^ = \/2mhh'i (a) Control ramps A(t). Red: original control 
ramp as derived from OCT. Black: control ramp after applying the electronics 
filtering. Blue: filtered, rescaled and shifted control ramp, (b) GPE momentum 
distribution, simulated without accounting for finite bandwidth. (c) GPE 
momentum distribution, simulation including finite bandwidth, rescaling of the 
control ramp by a factor of 1.6 and a time shift of 0.08 ms. 



3.3. Effect of finite control bandwidth 

The current in the modulation wire is driven by a custom-design low-noise current 
source, which is controlled from an arbitrary waveform generatoi[§j that outputs 
the excitation ramp. A slight smoothing of the control sequence is imposed by 
finite bandwidth of the electronics, which has to be accounted for when comparing 



experimental and numerical results (see Sec. 4.1). The measured transfer function 
modulus at a frequency v can be approximated by an exponential |A^(i^)| ~ 

^^/^co with cutoff frequency ~ 4.4 kHz. Furthermore, a frequency-dependent phase 
shift is imposed. Effectively, filtering causes a reduction of the driving amplitude near 
the resonant frequency vi ~ 1.8kHz by a factor |A^(^y)|~^ ~ 1.6, and a time delay 
on the order of 0.1ms (Fig. [sji). In Fig. |5]it is shown, that the filtering due to the 
electronics can be largely canceled by rescaling and shifting the control sequence by 
these factors. The difference in the outcome of the simulated momentum distribution is 



only small and largely given by a slightly enhanced collective oscillation (see Sec. 4.3). 



3.4' Measurement of momentum distributions 

At a time t after starting the excitation process, we suddenly switch off the trapping 
potential, implying that for t < 5 ms the excitation process is still incomplete. The 
fast transverse expansion of the cloud due to the tight waveguide confinement causes 
atom interactions to vanish rapidly, and the ensuing expansion can be considered 
ballistic. After texp = 46 ms of expansion, we take a fluorescence image [46j similar 
to that shown in Fig. [6]^a), fully integrating over the z-direction. In the images, three 
separate clouds can be observed along the longitudinal x-direction. The two side peaks 
emerge due to decay of the excited state which has been populated by our excitation 
protocol. Correlation properties of these (twin beams) and the dynamics of the decay 
process have been analyzed elsewhere [46, 13j. We separately integrate along x over 
the central and side peaks, respectively (blue dashed lines) to analyze the transverse 
state of each part of the system. The observed density distributions along y [Fig.[6]^b)] 
represent the momentary momentum distributions of the trapped cloud at time t, as 
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Figure 6. Experimental image analysis, (a) Typical experimental image 
data for optimal excitation and t = 5.5 ms, averaged over 12 shots. As the 
image is taken after 46 ms of expansion time, it predominantly reflects the initial 
momentum distribution. The scale bar corresponds to a spatial distance of 
187 iim, equivalent to the typical momentum of 5.5 pm"-*^ ko. (b) Transverse 
momentum distribution inferred from the image in panel (a). Upper (black) line: 
central peak, inside dashed lines in (a). Lower (red) line: emitted atoms, outside 
dashed lines in (a). Dashed (blue) line: Gaussian flt to emitted atom momentum. 



the initial transverse cloud size (of the order of ly ~ 250 nm) is negligible compared to 
that after expansion (far field). If we express momenta as wave numbers /Cy, a distance 
6y in the image hence corresponds to 6ky = a6y with a = m/httoi ^ 0.034 pm~^. 
Taking an experiment series where t is scanned, we can thus fully access the momentum 
distribution dynamics hyik^t) along the excitation direction, which we will typically 
depict as false-color plot, see e.g. Fig. [7|^a). Or main interest will be the dynamics 



of the central (source) cloud which is subject to the excitation; however, in Sec. 4.2 



also the transverse dynamics of the twin-beam peaks will be of some importance. For 
each of the experimental series shown in this paper, we averaged over typically 12 
realizations to suppress noise and allow for robust comparison to theory results. 



4. Results 



We will analyze the excitation dynamics in various complementary ways, motivated 
by the goal of developing an effective mapping of the many-body dynamics to a driven 
two-level system. In Sec. |4.1| we will start by comparing the results obtained in Sec. [2] 
for the time-dependent momentum distribution of the condensate wave function to 
experimental observations, varying a range of relevant parameters. While the excellent 
agreement ensures that the numerics used to obtain the optimized ramp are accurate, 
this result only gives limited insight into how the excitation process can be understood 
qualitatively. In Sec. |4.2| a more phenomenological analysis is performed directly on the 
experimental data, which will give hints about how to develop a two-level description. 
In Sec. |4.3| the GPE simulations are investigated in more detail, using a description 
based on Wigner quasi-probability functions, and displaced Fock states. It will become 
evident that all approaches lead to conceptually similar and quantitatively compatible 
interpretations, which can finally be unified to obtain a two-level interpretation as 
sought after initially. 
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Figure 7. Comparison of momentum distribution dynamics as obtained from 
experiment and theory for typical parameters, (a) Experiment. Each pixel column 
in the false-color plot corresponds to a distribution as shown in Fig.|6jb). (b) Id 
GPE numerics, including finite bandwidth effects (see text). 

4.1. Comparison of experiment and numerics 

Compared to other driven quantum systems, where optimal control techniques may be 
applicable, a rather unique advantage of cold atoms is the accessibility of the system 
response, enabled by the relatively large time and length scales and the abundance 
of powerful imaging techniques. Probing the performance of a control strategy such 
as that developed in Sec. [2] is not restricted to the final outcome, but the driven 
system can be monitored even while it is being driven, providing direct means to 
compare experiment and numerical simulations, or apply feedback schemes [47]. As 
explained in Sec. |3.4[ time-of-flight fluorescence images give us direct access to the 
momentum distribution of the condensate along its transverse axis. In Fig. [7|^a), a 
typical momentum distribution dynamics plot, as obtained from the experiment, is 
shown. 

Many-body effects Along the transverse directions, confinement is strong enough 
{hvi ^ ji) to make interaction-induced effects comparatively small. Still, to achieve 
the highest possible fidelity of the excitation, it is crucial to take into account the 
nonlinear term in Eq. ([T]) for optimization. In Fig. [sj the excitation dynamics is 
shown for a data set, where the atom number has been varied before starting the 
excitation sequence. The data is compared to the result of the GPE ([T]). 

It is observed, that effective excitation is achieved for a nonlinearity corresponding 
to an atom number N ~ 900, which is close to what has been used in the optimization. 
For all other atom numbers, stronger residual dynamics after the end of the sequence 
{t > 5 ms) is found, indicating decreased fidelity, as the desired state is stationary. 
While the GPE simulations reproduce the general tendencies found in the experiment, 
the agreement is not as good as e.g. for scaled excitations (see below). For the 
highest atom number, only rather poor qualitative agreement is reached, indicating 
insufficiency of a mean-field model such as GPE (necessitating e.g. a MCTDHB 
ansatz [121 128]) ^ind strong effects of the rapid decay of the excited state. 

Robustness against experiment inaccuracy In OCT, an aspect of high relevance is 
the sensitivity of the excitation dynamics to deviations of experimental parameters 
from the ones used for optimization. In our case, this predominantly applies to 
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Figure 8. Interaction effects on excitation. In each pair of plots (top: experiment, 
bottom: numerics), the typical experimental atom number is shown alongside the 
GPE nonlinearity term g (in Hz pm, ordinary frequency). The efficiency 77 is 
defined as described in the text. 
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Figure 9. Stability of the excitation sequence against inaccuracy of the trapping 
potential (numerical result). In each plot, the deviation of the potential terms 
Svy^Say are given (in units of Hz), as well as the efficiency 77 and the spurious 
excitation to higher states as defined in the text. 



parameters affecting the trapping potential. We consider small changes of the potential 
parameters Uy^ <jy. In the experiment such deviations arise from variation of the 
dressing parameters Aq, which, in turn, are caused by inaccuracy of the current 
in the RF antenna wires, and of the external offset field along x (defining the atomic 
Larmor frequency), respectively. 

Numerical results for a range of parameters are shown in Fig. [oj Panels (c-f) 
correspond to deviations caused by an offset field misalignment of ±2mG (b,c) and 
±7mG (e,f) leading to weaker (positive values) or stronger RF dressing, respectively. 
It is observed that any deviation leads to a decrease in excitation efficiency, which 
is defined here as time- averaged overlap with the desired wave function ^jj^^ rj = 
(I (?/^d|'0(^)) P)t>5ms- The similarly defined population of higher excited states ^ 
becomes strong at trap modifications with weaker dressing 5vy > and 5 ay < 0. This 
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effect can be expected, as the protection against excitation to higher states fades with 
decreasing anharmonicity, while the excursion of the trap relative to the typical length 
ly increases. In panel (d), on top of an offset field mismatch of +3mG, the current 
in the RF wire has been adjusted to cancel the effect on Vy. The weak mismatch in 
(jy and only leads to a slight reduction of efficiency. Consequently, optimizing the 
experimental parameters for a strong excitation (e.g. by minimizing non-stationarity 
at t > T = 5 ms) may lead to slightly shifted values, which however compensate. 
Using this method, a sensitivity of better than 1 mG (or an equivalent mismatch 
of the dressing current) can be reached, which is beyond what can be achieved by 
independent characterization of the trapping potential. 

Scaled excitations In Fig. [lO] the momentum distribution dynamics is shown for a 
data set with varying excitation efficiency, which will be the main subject of analysis 
in the remainder of this and the following section, as it covers a very broad range 
of control sequences. To achieve different efficiencies, the excitation ramp has been 
scaled in amplitude by factors s with respect to the optimal control result, resulting in 
strongly varying wave function dynamics. The approach of simple amplitude scaling 
has been chosen over using separately optimized ramps for different efficiencies, to 
allow for easier comparison due to the well-defined relation between the used control 
sequences. Furthermore, our analysis will show that the main spurious effect of this 
strategy are collective oscillations at reduced scalings. Comparison between GPE 
and experimental result (average over ~ 12 realizations) shows excellent agreement at 
early times. [[[] At later times, decay of the excited state into twin beams, which is 
not accounted for in theory, becomes significant (see bottom right panel) and for high 
values of s, agreement is reduced due to inelastic collisions with the twin beams which 
reside in a different transverse state. However, for weak excitation, even the shape of 
single "beating peaks" after the end of the excitation pulse is precisely captured by 
numerics. Along the /c-axis, the GPE result has been convolved with a Gaussian of 
^/{^hof) • 40 pm ^ 1.20 pm~^ rms width to account for finite imaging resolution and 
bulk position fluctuations. Apart from a small shift of the t-axis and a slight re-scaling 
of the /c-axisj5]the scaling factor s is the only free input parameter of the simulation. 

Having established the detection method, and verified that the outcome is 
consistent with the numerics on which the control optimization has been founded, 
we now proceed to a more qualitative analysis of the experimental result. 

4-2. Analysis of experimental momentum dynamics 

In this and the following section we will analyze the momentum distribution dynamics 
beyond a simple comparison to numerical results. The notion underlying the discussion 
will be that of a few-level system, comprised by the ground, first and occasionally 
second excited state of the confinement potential along the excitation direction, with 
the final goal to reduce the anharmonic oscillator to a closed two-level systemj^ This 

II Note that s has been defined including the necessary re-scahng due to finite electronics bandwidth 
(see Fig. Isl. 

^ The sniit in t is well below the experimental time resolution, and is very probably due to the 
inaccuracy of the filtering circuit characterization. The necessity for the re-scaling of k (of the order 
of 10%) might arise from interaction effects causing weak hydrodynamic effects in expansion j48] . 
The values of both adjustment are consistent among all sets shown. 

+ In the literature on quantum control strategies, this problem is occasionally discussed as that of 
leakage-suppression of a two- level system [49 50] 1191 [20] , 
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Figure 10. Results for scaled excitation ramps. Mean atom number is 770 for 
sets I, II, IV- VI, and 856 for sets III and VII. For each of the seven sub-sets, the 
upper image (red false-color) is the experimental result, normalized separately for 
each time step. The middle image (blue false-color) shows the numerical GPE 
result, including low-pass filtering and scaling by the factor s as given. The 
bottom image shows the deviation between experiment and theory, expressed as 
imbalance fiex — ^thJ the color scale for the imbalance is enhanced by a factor 
3. The bottom right inset shows the relative amount A^dec/-^ of atoms that have 
decayed from the excited state into twin atom pairs. 
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approach may seem inappropriate, as it relies on the superposition principle, which 
requires a linear equation of motion and is hence not applicable to a mean-field wave- 
function as described by the GPE. However, in our case the nonlinearity is weak 
compared to the oscillator energy, and so is the modification of the dynamics due to 
many-body effects (see Fig.|8|, suggesting that a description in terms of single-particle 
states may still provide significant insight. 



Center- of -mass dynamics As the simplest possible observable derivable from the 
momentum dynamics, we start by analyzing the transverse center-of-mass of the 
experimental images, corresponding to the momentum expectation value K{t) = 



(/Cy(t)), see black lines in Fig. 11 (left panels). In the power spectra of K{t) (center 
panels), two strong peaks are observable near the first two transverse level spacings at 
frequencies ui = 1.84 kHz and U2 = 1.99 kHz, and a weak third at z/3 ^ 2.10 kHz, 
defined analogously. Assuming a single-particle level picture, these peaks can be 
interpreted as beating frequencies between populations of the first three levels of the 
oscillator, where mean-field effects are causing frequency shifts, as described below. 
Consequently, the magnitude of oscillations is the strongest for intermediate excitation 
efficiencies (sets II-IV), where the levels are populated most evenly, maximizing the 
beating contrast (see below). 

A crucial observation is, that also the transverse profiles of the twin-beam peaks, 
which are separated in the images longitudinally (see Fig.[6|, exhibit strong oscillations 
of Kt{t) = {ky\t)). Meanwhile, they fully maintain their Gaussian shape [Fig. [oj^b)]. 
In Fig. [11] oscillations of the relative center-of-mass Kj.{t) = K{t) — K^{t) (left), and 
their power spectrum (center) f{v) = \T[Kr{t)]{i')\'^ , are shown as blue lines. It is 
observed that, while the oscillations are similarly strong as in a fixed frame, all peaks 
in the power spectrum, except that near ui are suppressed. This suggests, that in a 
reference frame co-oscillating with Kt {t) , the dynamics can be understood in terms of 
two transverse levels, motivating an approach of decomposition into a quasi-classical 
oscillation, and "internal" dynamics, which remain unaffected by the bulk oscillationj^ 
This interpretation is consistent with our understanding of the decay process [M], 
where the transverse state of the twin beams, a ground state displaced by Kt{t), 
defines the appropriate ground state for the internal dynamics. In Sec. |4.3[ a more 
rigorous formalism for the co-oscillating frame will be given, and its position will be 
independently derived from numerical results. 

In the right column of Fig. [TT] spectra are shown which are derived from the 
oscillations at times t > 5 ms only, i.e., where no driving occurs anymore. Hence, 
they provide a characterization of the final state that is reached after the excitation. 
Qualitatively, the same features are observed as in the full time spectra, however, 
peaks at 1^2 are smaller, which is consistent with theory, as will be shown below. Also, 
in the relative center-of-mass spectrum, the observation of a single-peak structure, 
with a minimal amplitude for the most efficient excitation is even more evident. 



In Fig. 12 ^a), the integrated power of the oscillations P cx J f{u)du^ measuring 
the stationarity of the final state, is shown as a function of the numerically obtained 
excitation efficiency r] (see previous section). Apart from the strongest driving, where 
higher states may become excited more easily, P shows fair agreement with a curve 

* This decomposition is exactly valid for harmonically confined many-body systems [171 118] . 
Obviously, this does not hold for an anharmonic oscillator, which is exactly why our excitation to a 
non-classical state by displacement can work at all (see above). Being aware of the inconsistency, we 
still apply the decomposition approach to qualitatively understand the dynamics. 
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Figure 11. Momentum space center-of-mass dynamics for data set as shown in 
Fig. |10| Left column: center-of-mass momentum of the source cloud with respect 
to a fixed frame (K (t), black) and relative to the twin-beam center-of-mass (Kr(t), 
blue). (See Fig. |14| for the twin-beam center-of-mass.) In the background, the full 
dynamics is shown (see Fig. |10| )- Middle column: corresponding power spectra 
/(z/), taken over the entire time span shown. Right column: spectra, taken over 
a time span starting from t > T = 5 ms, i.e., after the end of the excitation. Grey 
dashed lines in the background indicate the harmonic frequency U]^, and the first 
three level spacings, as defined in the previous section. All spectra are in arbitrary 
units, but normalized identically for each of the columns. 



Vibrational state inversion of a BEC 



19 




Figure 12. Analysis of post-excitation beating spectra shown in the right column 
of Fig. [ll] (a) Integrated power of oscillations P. The experimental points have 
been scaled along the y-axis for best fit to 77(1 — 77) (red line). 77 has been derived 
as described in the previous section. The shading of each point indicates the 
corresponding scaling s (white is highest), (b) Peak position (black, left axes) 
and cosine of averaged phase (green, right axes). Red and blue lines correspond 
to the single-particle level spacing ui , and the mean-field-shifted level spacing u'^ , 
respectively. 



given by 7^(1 — 7^), which is the squared amphtude of the interference term in the 
momentum-space density of a two-level system with momentum- space wave functions 



= (1- 



^)IV^o(MI'- 



hv\Mky)\ 



(10) 



2y/r]{l-r]m*o{ky)Mky)] cos(2^^;t). 



The positions of the beating peak (obtained from a Gaussian fit) are shown in 
Fig, p^b). For high efficiency, the frequency is shifted downwards from the oscillator 
level spacing ui (red line). This is explained by the mean field term in the GPE ([T]). 
For the boundary case of near-unity efficiency, the shift can be calculated rather easily. 
As the ground state population is negligible, it does not contribute to the interaction 
energy, and the chemical potential fie for an atom in the excited state is given by 
the second eigenvalue of the time-independent GPE. The according wave function 
?/^d (i.e., the desired state in the optimization process) can now be used to calculate 
the chemical potential of a single atom in the ground state ipQ^ using a Schrodinger 
equation with effective potential arising from the mean field of the excited state: 



/ieV^o(^) = 



2m dy'^ 



VeAy)+2g\My)f 



V'o(2/)- 



(11) 



The beating frequency is now given by the difference in chemical potential. Instead of 
the oscillator level spacing ui ^ 1.831 kHz, we obtain u[ = (/ig— M^)/(27rfi) ~ 1.724 kHz 
(blue line). Given the uncertainty in the input parameters of the calculation (such 
as the assumption of an equilibrium Thomas-Fermi shape longitudinally), this value 
agrees well with the experimentally obtained one for maximum efficiency (set IV), 
= 1.709(4) kHz. 

Finally, we can have a look at the phase of the (relative) center-of-mass oscillation. 
When comparing the value of Ky'^ (t) for different scalings at a fixed time in Fig. [ll] 
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it is apparent, that the phase inverts at the point of maximum efficiency. We take the 
averaged phase from the Fourier transform result, weighted by the Lorentzian fit of the 



peak, and obtain the curve shown in Fig. 12 ^b, right axes). The inversion is reminiscent 



of a two-level system subject to a Rabi driving, where the phase of precession inverts 
after passing the pole of the Bloch sphere at a pulse area larger than tt. As will be 
shown in Sec. |4.3[ the excitation process can be understood analogously. 

4.3. Interpretation of numerical result: two-level driving model 

To understand the physical mechanism governing the optimal excitation protocol, in 
the following we analyze the time evolution of the condensate wave function in the 
Wigner representation [16j : 

W{y,k,t) = le-"^^^{y+'-,t)r{y-^,t)ds, (12) 

which provides a mixed position-momentum distribution. The Wigner function has 
many appealing features reminiscent of a classical distribution function. Integration 
over all momenta k gives the spatial probability distribution Likewise, 
integration over y gives the momentum probability distribution. The Wigner function 
of the condensate ground state. Fig. [Tsj^f), approximately corresponds to the ground 
state of the harmonic oscillator, with equal uncertainty in position and momentum. 
In the figure the distribution is slightly elongated along y due to the nonlinear atom- 



atom interactions. The desired state of the control. Fig. 13 ^b), corresponds to the first 
excited state of the GPE in the anharmonic trap. It has positive and negative values 
(giving a node ai y = upon integration over all momenta), and thus differs from a 
genuine classical distribution function. 

Panel (a) of the figure reports the time evolution of the Wigner function. We 
plot the iso-surfaces at ±0.35 times the maximum value of the Wigner function. At 
times later than 3 ms we have added transparency to the iso-surface for positive 
values to show the appearance of the negative part of the Wigner function, associated 
with a non-classically excited state. The solid line shows the time variation of the 
spatial minimum of the confinement potential. Initially, this displacement brings the 
condensate into collective oscillations, whose frequency is determined by the harmonic 
part of the confinement potential. For a large enough displacement, the condensate 
wave function is brought into the region where the anharmonicity of the confinement 
is sufficiently large to modify the internal structure of the wave function (and not just 
its displacement). One observes that in addition to the center-of-mass oscillations in 
this regime the transfer from the ground to the first excited state occurs. Finally, at 
the terminal time T = 5 ms of the control process the minimum of the confinement 
potential is shifted to bring the condensate to a complete halt. 

We next suggest a procedure to approximately map the excitation dynamics onto 
a genuine two- level description of ground and excited condensate states. As in Sec. |4.2[ 
the main idea is to separate the wave function dynamics into (i) a collective, quasi- 
classical oscillation, which is needed to bring the condensate into the anharmonic part 
of the trap, and (ii) an internal conversion between the ground and first excited state, 
defined in a co-moving frame. The latter conversion is governed by the anharmonic 
part of the trap, as explained above. 

We define wave functions (/>g(y) and (l)e{y) as single-particle eigenfunctions of 
the harmonic part of the trap potential only, i.e. Eq. ^ with ay and set 
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Figure 13. Time evolution of Wigner function. Panel (a) reports the time 
evolution of the Wigner function, and panels (b-f) show snapshots at selected 
times. In (a) we show the iso-surfaces at zbO.35 times the maximum value of the 
Wigner function, with transparency added to the iso-surface at time above 3 ms 
to highlight the appearance of the negative Wigner function part, associated with 
the first excited state. For discussion see text. 
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to zero. Also, any modifications due to the nonlinear atom-atom interactions are 
neglected. This simplification allows us to analyze the dynamics in terms of displaced 
Fock states [51], that capture well the notion of the separation approach. Let 
D[a{t)\ = exp[a{t)a^ — a{t)*a] denote the displacement operator of the harmonic 
oscillator [16j, where a{t) = [l~^YQ{t) + ilyKo{t)]/y/2 determines the position and 
momentum of the displacement at time and d denotes the annihilation operator. 
For a given displacement we can compute the overlap between the displaced 

ground and excited states with the condensate wave function according to 



(13) 



Determining the value a{t) which gives the largest overlap at time t allows us the 
aforementioned decompositions into (i) center-of-mass coordinates Yo{t) and Ko{t)^ 
and (ii) probability amplitudes {D[a{t)](j)g\ip{t)) and {D[a{t)](j)e\ip{t)) for the ground 
and excited state within the displaced frame. In all cases we find an overlap x{t) 
well above 90%, which thus justifies the wave function decomposition. In Fig. [l4j 
the obtained values for Ko{t) are shown as red lines, and compared to experimentally 
obtained values, as described below. In Fig. [Ts] we compare x(^) and the obtained 
excited population r]'{t) = [{D[a{t)](j)g\ilj{t))fto results from direct projection of 
il^{y^t) on the oscillator states (/)(?/) (defined in the co-moving frame of the excitation 
motion). The direct projection leads to strong transient population of higher excited 
states, and a sudden jump near the end of the excitation [Fig. [Tsj^a)], where they 
are depopulated again. This is reminiscent of the fixed-frame center-of-mass spectra 
(black lines in Fig. [TT] ), where a peak near V2 is present when regarding the entire 
sequence (center column), but mostly vanishes after t = T (right column). In 
contrast, the two- level approximation in the system displaced by a{t) yields a smooth 
transition [Fig. [T5|b)], consistent with the continuous appearance of negative values 
of the Wigner function (Fig. 13). In the momentum dynamics derived from the two- 
level model in a similar manner to Eq. (10) as shown in Fig. flSld), a continuous 



transfer to the excited state is observed, with strong beating at intermediate excited 
population. Again, this is consistent with the experimental relative center-of-mass 



spectra (blue lines in Fig. 11), where only a single peak near vi persists, even during 
the excitation. Similar to a Rabi pulse with area larger than tt, the excited population 
rj{t) is decreasing towards t = T for scaling parameters 5 > 1. 

As laid out in Sec. |4.2[ the appropriate ground state for the internal conversion 
dynamics can also be determined in the experiment from the center-of-mass position 
Kt{t) of the twin-atom beams which the excited state is decaying into continuously. 
For times t, where the decayed fraction becomes perceivable, we can compare the 



experimentally found Kt{y) to Ko{y) as in Fig. 14, and find good agreement without 
any free parameter over a large range of settings. Together with the absence of decay 
products from higher excited states in the experiment, this result confirms the validity 
of the decomposition approach. In Ref. [M] it has been shown, that the obtained 
populations of the excited state lead to an accurate quantitative description of the 
ensuing decay process. 



5. Conclusion 



In conclusion, we have presented successful application of optimal control theory to the 
problem of preparing a non-classical, strongly out-of-equilibrium motional state of a 
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Figure 14. R efer ence frame for two-level model. Underlying data are the same as 
shown in Figs. p^ and |ll[ Red lines are the momentum-space displacement Ko(t) 
of the two-mode basis states, as obtained from applying Eq. \13\ to the GPE 
result. Black points indicate the experimentally found center-of-mass position of 
twin beams that have decayed from the excited state Kt, defining the reference 
frame for the emission process (see Sec. |4.2| ). Similar to the momentum space 
dynamics as shown in Fig. |10| agreement reduces at later times, where decay into 
twin beams becomes strong. Data set I has been omitted due to the emission of 
twin beams being insufficient to determine Kt. 
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Figure 15. State populations during the excitation process, (a) Populations 
?7(t) of excited states of the full anharmonic potential [Eq. ([7])] as arising from 
direct projection of the GPE result for data set V (s = 1). The solid line 
indicates the population of the first excited state, dotted lines represent the 
ground (black) and first and second excited (red, blue) states, (b) Corresponding 
momentum dynamics (identical to Fig. |10| -V). (c) Population of first excited state 
in co-oscillating frame within the two-mode model r]' {€). Solid line: data set V, 
corresponding to solid line in (a). Dashed lines: sets I (black, s = 0.27), III (red, 
s = 0.63), and VII (blue, s = 1.43). The dash-dotted line indicates the total 
overlap of the two- level model with the GPE result x(0 [^^^ (|13|>], which 
exceeds a value of 0.95 at all times t. (d) Momentum dynamics arising from time- 
dependent superposition of (f)o,(j)i in the co-oscillating frame, data set V. Note 
the strong beating at intermediate times/excited fractions. 
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Bose-Einstein condensate, realizing population inversion with near-unity fidelity. The 
obtained condensate wave function corresponds to the first excited eigenstate of the 
Gross-Pitaevskii equation, closely resembling the first odd Fock state of a harmonic 
oscillator. To manipulate the external state of the Bose-Einstein condensate, we 
used precisely controlled motion of an anharmonic trap potential along the optimized 
trajectory. Experimental and numerical results on the momentum distribution 
dynamics during and after the excitation sequence show excellent agreement over a 
large range of parameters, including tuning of many-body effects. Moreover, a model 
of the excitation dynamics based on decomposition into a quasi-classical oscillation 
and the actual state transfer has been developed, and shown to be consistent with 
various observations made in both experiment and theory. Using this approach, we 
were able to deduce an approximate two- level description of the excitation process. 

A first application of the vibrational state inversion, using the condensate as a 
gain medium for matter wave amplification, has been demonstrated in Refs. [T3l [T4]. 
However, optimal control in condensates is not restricted to high-fidelity preparation 
of a desired wave function, and more general pulses that e.g. acting on non- 
stationary initial states in a phase-sensitive manner can be implemented [52] , 
State preparation beyond a mean-field description has been proposed, including 
entanglement generation [53) [54j. number-squeezed states [12j, or cooling [55j, which 
should be realizable in a similar fashion. More generally, our results highlight the 
potential of experiments with Bose-Einstein condensates as a test-bed for a large 
range of quantum control problems, as known from NMR spectroscopy [56l |57l [58] 
solid-state, [591 [60l Ell [20l |50l [19] , atomic, [62j, or molecular physics p3l 164} [65] . 
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